Local Density Fluctuations, Hyperuniformity, and Order Met 



Salvatore Torquato 
Department of Chemistry and Princeton Materials Institute, 
Princeton University, Princeton, NJ 08544 

Frank H. Stillinger 

Department of Chemistry, Princeton University, Princeton, NJ 08544 

(Dated: February 2, 2008) 



Abstract 

Questions concerning the properties and quantification of density fluctuations in point patterns 
continues to provide many theoretical chahenges. The purpose of this paper is to characterize 
certain fundamental aspects of local density fluctuations associated with general point patterns in 
any space dimension d. Our specific objectives are to study the variance in the number of points 
contained within a regularly-shaped window J7 of arbitrary size, and to further illuminate our un- 
derstanding of hyperuniform systems, i.e., point patterns that do not possess infinite- wavelength 
fluctuations. For large windows, hyperuniform systems are characterized by a local variance that 
grows only as the surface area (rather than the volume) of the window. We derive two formula- 
tions for the number variance: (i) an ensemble-average formulation, which is valid for statistically 
homogeneous systems, and (ii) a volume-average formulation, applicable to a single realization of a 
general point pattern in the large-system limit. The ensemble-average formulation (which includes 
both real-space and Fourier representations) enables us to show that a homogeneous point pattern 
in a hyperuniform state is at a "critical-point" of a type with appropriate scaling laws and crit- 
ical exponents, but one in which the direct correlation function (rather than the pair correlation 
function) is long-ranged. We also prove that the nonnegativity of the local number variance does 
not add a new realizability condition on the pair correlation. The volume-average formulation is 
superior for certain computational purposes, including optimization studies in which it is desired to 
find the particular point pattern with an extremal or targeted value of the variance. We prove that 
the simple periodic linear array yields the global minimum value of the average variance among all 
infinite one-dimensional hyperuniform patterns. We also evaluate the variance for common infinite 
periodic lattices as well as certain nonperiodic point patterns in one, two, and three dimensions 
for spherical windows, enabling us to rank-order the spatial patterns. Our results suggest that the 
local variance may serve as a useful order metric for general point patterns. Contrary to the con- 
jecture that the lattices associated with the densest packing of congruent spheres have the smallest 
variance regardless of the space dimension, we show that for d = 3, the body-centered cubic lattice 
has a smaller variance than the face-centered cubic lattice. Finally, for certain hyperuniform disor- 
dered point patterns, we evaluate the direct correlation function, structure factor, and associated 
critical exponents exactly. 
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I. INTRODUCTION 



The characterization of density fluctuations in many-particle systems is a problem of great 
fundamental interest in the physical and biological sciences. In the context of liquids, it is 
well known that long- wavelength density fluctuations contain crucial thermodynamic and 
structural information about the system |l| . The measurement of galaxy density fluctuations 
is one of the most powerful ways to quantify and study the large-scale structure of the 
universe . Knowledge of density fluctuations in vibrated granular media has been used 
to probe the structure and collective motions of the grains Ij]. Recently, the distribution 
of density fluctuations has been employed to reveal the fractal nature of structures within 
living cells ^. 

Clearly, density fluctuations that occur on some arbitrary local length scale 
llol | provide considerably more information about the system than only long-wavelength 
fluctuations. Our main interest in this paper is to characterize certain fundamental aspects 
of local density fluctuations associated with general point patterns in any space dimension 
d. The point patterns may be thought as arising from the coordinates of the particles in a 
many-particle system, such as the molecules of a liquid, glass, quasicrystal, or crystal, stars 
of a galaxy, grains of a granular packing, particles of a colloidal dispersion, or trees in a 
forest. 

Consider an arbitrary point pattern in (i-dimensional Euclidean space 3?*^. Let Q represent 
a regular domain (window) in S?*^ and xq denote a conflgurational coordinate that specifles 
the centroid of the window Q. The window will always have a flxed orientation. There 
is a variety of interesting questions that one could ask concerning the number of points 
contained within Q. For example, how many points Nq are contained in Q at some flxed 
coordinate xq? This question is a deterministic one if the point pattern is regular and may 
be a statistical one if the point pattern is irregular (see Fig. Q). How does the number of 
points contained within some initially chosen Q at flxed coordinate xq vary as the size of 
Q is uniformly increased? How do the number of points within a flxed Q fluctuate as xq is 
varied? 

For a Poisson point pattern, the statistics of the number of points contained within a 
regular domain are known exactly. For example, the number variance is given by 

{N'n)-{Nny = {Nn), (1) 



FIG. 1: Schematics indicating a regular domain or window 17 and its centroid xq for two different 
point patterns. Left panel: A periodic point pattern. Right panel: An irregular point pattern. We 
will show that the statistics of the points contained within for these two types of patterns are 
fundamentally different from one another. 

where angular brackets denote an ensemble average. Letting be a (i- dimensional sphere 
of radius R and noting that (Nq) is proportional to R^, leads to the result that the number 
variance grows as the sphere volume, i.e., 

(iV2) - {Nnf oc R'. (2) 

This result is not limited to Poisson point patterns. Indeed, a large class of correlated 
irregular point patterns obeys the variance formula as we will discuss in Section II. 

Can the variance grow more slowly than the volume of the domain or window? One can 
show that for any statistically homogeneous and isotropic point pattern, the variance cannot 
grow more slowly than the surface area of the domain, whether it is spherical or some other 



strictly convex shape 
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Q| . Thus, it is natural to ask the following question: For what 
class of point pattern does the variance grow as the surface area? For a spherical domain, 
we want to identify the point patterns that obey the variance relation for large R 

(iV^) - (Nn)' ~ R"-'. (3) 

We will refer to such point patterns as "hyperuniform" systems because, as we will see, 
such systems do not possess infinite-wavelength fluctuations. (This is to be contrasted with 
"hyposurficial" systems, whose "surface" fluctuations vanish identically.) Additionally, it 



is of great interest to identify the particular point pattern that minimizes the amphtude 
(coefficient) of the fluctuations that obey or achieves a targeted value of this coefficient. 

Clearly, points arranged on a regular (periodic) lattice are hyperuniform. More generally, 
it is desired to know how the number of lattice points N{R) contained within a spherical 
window of radius R varies as function of R when the sphere is centered at xg. For simplicity, 
let us consider this question in two dimensions for points arranged on the square lattice and 
let the center of the circular window of radius R be positioned at a point (ai, 02) in the unit 
square. The answer to this query amounts to finding all of the integer solutions of 



(ni - aif + {n2 - aa)^ < 



(4) 



m 



IJ]. This problem is directly related to the 



a problem of interest in number theory 
determination of the number of energy levels less than some fixed energy in integrable 
quantum systems ^. It is clear that N{R) asymptotically approaches the window area 
ttR^ and unit density, for large R. The apparent "random" nature of N{R) is beautifully 
illustrated in Figure which shows how the function N{R) — ttR^ grows with R. 
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FIG. 2: The function N{R) — -kB? versus R for the unit-spacing square lattice, using a circular 
window of radius R centered on a lattice point. 



It is considerably more challenging to identify non-periodic point patterns, such as dis- 
ordered and quasiperiodic ones, that are hyperuniform. The mathematical conditions that 
statistically homogeneous hyperuniform systems must obey (derived in Section II) are a 
necessary starting point in identifying such hyperuniform point patterns. These conditions, 
which include the counterintuitive property of a long-ranged "direct" correlation function. 
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are determined from a general formula for the number variance of such systems, which is ob- 
tained in Section II. The fact that the direct correlation function of a hyperuniform pattern 
is long-ranged is reminiscent of the behavior of the pair correlation function of a thermal 
system near its critical point. Indeed, we show that a statistically homogeneous point pat- 
tern in a hyperuniform state is at a "critical-point" of a type with appropriate scaling laws 
and critical exponents. By deriving a Fourier representation of the local variance, it is also 
shown that the nonnegativity of the variance does not add a new realizability condition on 
the pair correlation function beyond the known ones. 

To date, only a few statistically homogeneous and isotropic patterns have been rigorously 
shown to be hyperuniform. One of the aims of this paper will be to identify other such 
hyperuniform examples, and to describe a procedure to find them systematically. This 
requires a formulation for the local variance that can be applied to a single realization of 
any pattern, which is accomplished in Section III. In Section IV we prove that the simple 
periodic linear array yields the global minimum value of the average variance among all 
infinite one-dimensional hyperuniform patterns. Interestingly, we also show that the variance 
for large spherical windows enables us to rank-order common regular lattice and certain 
disordered point patterns in one, two, and three dimensions (sec Sections IV and V). Our 
results suggest that the local variance may provide a useful order metric for general point 
patterns (see Section VI). Contrary to the conjecture that the Bravais lattice associated 
with the densest packing of congruent spheres has the smallest variance regardless of the 
space dimension, we show that for d — 3, the body-centered cubic lattice has a smaller 
variance than the face-centered cubic lattice. In Section V, we evaluate the direct correlation 
function, structure factor, and associated critical exponents exactly for certain hyperuniform 
disordered point patterns. Three appendices provide analytical formulas for key geometrical 
quantities required for the theory, an evaluation of the variance for hard rods in equilibrium 
for large windows, and a discussion of a certain property of hyposurficial point patterns. 

II. LOCAL VARIANCE FORMULA FOR REALIZATIONS OF STATISTICALLY 
HOMOGENEOUS SYSTEMS 

A general expression for the local number variance for realizations of statistically homo- 
geneous point patterns in d dimensions is derived. This is necessarily an ensemble-average 
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formulation. We obtain both a real-space and Fourier representation of the variance. Prom 
these results, we obtain formulas for asymptotically large windows. We show that a hy- 
peruniform point pattern is at a type of "critical-point" with appropriate scaling laws and 
critical exponents, but one in which the direct correlation function is long-ranged. 



A. Preliminaries 

Consider N points with configuration = ri, r2, . . . , rjv in a volume V. The local 
number density at position x is given by 

N 

n(x) = ^5(x-r,), (5) 
1=1 

where 5(x) is the Dirac delta function. The point pattern is statistically characterized by 
the specific probability density function P/v(r^), where P]s[{r^)dr^ gives the probability 
of finding point 1 in volume element dri about ri, point 2 in volume element dr2 about 
r2, . . ., point in volume element drjy about rjv- Thus, Pjvli"^) normalizes to unity 
and dr^ = dri, dv2, . . . , dr^v represents the A^d-dimensional volume element. The ensemble 
average of any function f{r^) that depends on the configuration of points is given by 

^^^'""^^ = IvIv---Iv (6) 
Because complete statistical information is usually not available, it is convenient to in- 
troduce the reduced generic density function p„(r") {n < N), defined as 

where rfr^"" = (ir„_|_i(ir„_|_2 ■ ■ ■ dr]\[. In words, p„(r")(ir" is proportional to the probability of 
finding any n particles (n < N) with configuration r" in volume element dr"-. In light of its 
probabilistic nature, it is clear that Pn(r") is a nonnegative quantity, i.e., Pn(r"') > 0, Vr". 

For statistically homogeneous media, Pni'"^"') is translationally invariant and hence de- 
pends only on the relative displacements, say with respect to ri: 

p„(r") = p„(ri2,ri3, . . . ,ri„), (8) 

where r^j = rj — r^. In particular, the one-particle function pi is just equal to the constant 
number density of particles p, i.e., 

N 

pi(rO=p^^hm^-. (9) 
7 



The limit indicated in (jHj) is referred to as the thermodynamic limit. Since our interest in 
this section is in statistically homogeneous point patterns, we now take the thermodynamic 
limit. It is convenient to define the so-called n-particle correlation function 

9n{rn = (10) 

In systems without long-range order and in which the particles are mutually far from one 
another (i.e., r^j = \rij\ — oo, 1 < i < j < N), Pnij'^) p" and we have from (|Tnjl that 
(7„(r") 1. Thus, the deviation of Qn from unity provides a measure of the degree of spatial 
correlation between the particles, with unity corresponding to no spatial correlation. 
The important two-particle quantity 

/ N P2(ri2) 

^2(ri2) = — (11) 

P 

is usually referred to as the pair correlation function. The total correlation function h{ri2) 
is defined as 

h{ru) = ^72(ri2) - 1, (12) 

and thus is a function that is zero when there are no spatial correlations in the system. When 
the system is both statistically homogeneous and isotropic, the pair correlation function 
depends on the radial distance ri2 only, i.e., 

5'2(ri2) = 5'2(ri2), (13) 

and is referred to as the radial distribution function. From fjll|) . we see that psi{r)g2{r)dr 
is proportional to the conditional probability of finding a particle center in a spherical shell 
of volume si{r)dr, given that there is another at the origin. Here Si(r) is the surface area 
of a d-dimensional sphere of radius r, which is given by 

where r(x) is the gamma function. Hence, for a finite system, integrating (A^ — I)g2{r)/V 
over the volume yields — 1, i.e., all the particles except the one at the origin. 

Observe that the structure factor 5'(k) is related to the Fourier transform of /i(r), denoted 
by /i(k), via the expression 

S(k) = l + p/i(k). (15) 
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The Fourier transform of some absolutely integrable function /(r) in d dimensions is given 

by 

/» = J f{r)e-'^-Ur, (16) 

and the associated inverse operation is defined by 

/(r)=(^//(k)e^''-'^k, (17) 

where k is the wave vector. It is well known that the structure factor is proportional to 
the scattered intensity of radiation from a system of points and thus is obtainable from a 
scattering experiment. An important property of the structure factor is that it must be 
nonnegative for all k, i.e., 

S{k) > Vk. (18) 



B. General Variance Formulas 



Let R symbolize the parameters that characterize the geometry of the window il. For 
example, in the case of an ellipsoidal window, R would represent the semi-axes of the 
ellipsoid. Let us introduce the window indicator function 

for a window with a configurational coordinate Xq. The number of points Nq within the 
window at Xq, which we henceforth denote by A'"(xo; R), is given by 



A^(xo;R) = / n(x)w(x — xq; R)(ix 

N 

= Yl 5(x-rj)w(x-xo;R)(ix 

N 

= ^w;(r,-xo;R). (20) 
1=1 

Therefore, the average number of points contained within the window in a realization of the 
ensemble is 

N 

(iV(R)) = / ^w(r,-Xo;R)P^(Odr^ 

•'^ i=l 

= / pi(ri)w(ri - xo;R)ciri 
— pi voir-. R)dr 

= P^i(R-), (21) 



where vi (R) is the volume of a window with geometric parameters R. Note that translational 
invariance of the point pattern, invoked in the third hne of relation (pT|. renders the average 
(A^(R)) independent of the window coordinate xq. 

Similarly, ensemble averaging the square of (j!2(Jj) and using relation (j2H) gives the local 
number variance as 

(iV2(R)) - (iV(R))2 = f pi(ri)^(ri-xo;R)rfri 

J J [p2(ri,r2) -pi(ri)pi(r2)]w(ri - Xo;R)w(r2 - xo; R)rfrirfr2 



+ 



(iV(R)) 



p / h(r)a(r; iVjdr 



where h{r) is the total correlation function defined by ()12j) . 



(22) 



and 

<(r; R) = / w(ri - xq; R)w(r2 - xq; R)rfxo (24) 

is the intersection volume of two windows (with the same orientations) whose centroids 
are separated by the displacement vector r = ri — r2 Appendix A provides explicit 

analytical formulas for the intersection volume for spherical windows in arbitrary dimension 
d. As before, statistical homogeneity, invoked in the second line of (|^^. renders the variance 
independent of xq. 

Remarks: 

n 

1. Formula (j22j) was previously derived by Landau and Lifschitz 16], although they did not 
explicitly indicate the scaled intersection volume function a(r;R). Martin and Yalcin 17 1 
derived the analogous formula for charge fluctuations in classical Coulombic systems. 

2. The local variance formula (1^^ is closely related to one associated with the local volume 

nn 

fraction fluctuations in two-phase random heterogeneous materials [15|, 118|| . Both formulas 
involve the scaled intersection volume function a{r; R). The essential difference is that the 
variance for local volume fraction fluctuations involves a different correlation function from 
h{r), namely, the probability of finding two points, separated by a displacement r, both in 
the same phase. 

3. The existence of the integral in ()22|) requires that the product /i(r)a(r; R) be integrable. 
For finite size windows, this will be the case for bounded h{r) because a(r; R) is zero beyond 
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a finite distance. For infinitely large windows, a(r; R) = 1, and integrability requires that 
h{r) decays faster than |r|~'^"'"^ for some e > 0. For systems in thermal equilibrium, this will 
be the case for pure phases away from critical points. The structure factor S(k) [defined by 
(jl5|) ] at k = diverges as a thermal critical point is approached, implying that h{r) becomes 
long-ranged, i.e., decays slower than |r|~'^ [19j. 

An outstanding question in statistical physics is: What are the existence conditions for 
a valid (i.e., physically realizable) total correlation function h{r) [2^ of a point process at 
fixed density p? The generalization of the Wiener-Khinchtine theorem for multidimensional 

n 

spatial stochastic processes ^21] states a necessary and sufficient condition for the existence 
of an autocovariance function of a general stochastically continuous homogeneous process is 
that it has a spectral (Fourier-Stieltjes) representation with a nonnegative bounded measure. 
If the autocovariance is absolutely integrable, this implies that its Fourier transform must be 
nonnegative. The total correlation function h{r) is the nontrivial part of the autocovariance 
function for a point process, i.e., it excludes the delta function at the origin. The fact that 
h{r) comes from a statistically homogeneous point process, however, would further restrict 
the existence conditions on h{r) beyond the Wiener-Khinchtine condition, which amounts 
to the nonnegativity of the structure factor. Obviously, besides the condition that 5'(k) > 0, 
we have the pointwise condition h{r) > — 1 for all r. The determination of other realizability 
conditions on h{r) is a open question j2l|. 

Thus, it is interesting to inquire whether the nonnegativity of the local number variance, 
given by formula p2|l . is a new condition on h{r) beyond the nonnegativity of the structure 
'actor 5'(k). As we now prove, the answer is no. By Parseval's theorem for Fourier transforms 



22| |. we can rewrite the general variance formula (|22j) for an arbitrarily shaped (regular) 



window as 

{N\R)) - (iV(R))2 = {N{R)) 

where 



^+ (2^/ W«(k;R)t/k 



(25) 



«(k; R) = > (26) 

is the Fourier transform of the scaled intersection volume function ()23j) and w(k; R) is the 
Fourier transform of the window indicator function (fT!H) . Again, by Parseval's theorem 

^ ^6(k;R)rfk= /u;2(r)rfr = 1. (27) 



(27r)^y ' ' ' t;i(R) 
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Finally, utilizing the definition psp of the structure factor, we arrive at the Fourier repre- 
sentation of the number variance: 



(N'iR)) - (iV(R))^ = (iV(R)) 



s(k)a(k- R)dk 



(28) 



Interestingly, we see that the variance formula can be rewritten in terms of the structure 
factor and the nonnegative function a(k; R), the Fourier transform of the scaled intersection 
volume function a(r; R): a purely geometric quantity. Since the latter is independent of the 
correlation function h{r), we conclude that the nonnegativity of the number variance does 
not introduce a new realizability condition on h{r). 

Remarks: 

1. Given the Fourier representation formula ()28p . it is simple to prove that the local number 
variance is strictly positive for any fi(R) > 0. Both the functions Q;(k;R) and S'(k) are 
nonnegative. Therefore, because the nonnegative integrand of formula (j28|) cannot be zero 
for all k, it immediately follows that the local variance is strictly positive for any statistically 
homogeneous point pattern whenever fi(R) > 0, i.e., 

{N^(R)) - {N(R))^ > 0. (29) 



2. Let the window grow infinitely large in a self-similar (i.e., shape- and orientation- 
preserving) fashion. In this limit, which we will denote simply by fi(R) — > oo, the function 
Q;(k; R) appearing in ()28|) tends to (27r)'^(5(k), where 5(k) is a d-dimensional Dirac delta 
function, and therefore dividing the variance (j^Hj) by (A^(R)) yields 

Observe also that the form of the scaled variance ()30p for infinitely large windows (or infinite- 
wavelength limit) is identical to that for equilibrium "open" systems, i.e., grand canonical 
ensemble, in the infinite-system limit. It is well known that the variance in the latter instance 
is related to thermodynamic compressibilities or susceptibilities The important distinc- 
tion is that result ()3()|1 is derived by considering window fluctuations in an infinite "closed" 
possibly nonequilibrium system. When the point pattern comes from a statistically homo- 
geneous equilibrium ensemble, one can interpret the fluctuations as arising from differences 
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in the point patterns in the ensemble members for a fixed window position or, equivalently, 
from moving the asymptotically large window from point to point in a single system. The 
latter scenario can be viewed as corresponding to density fiuctuations associated with an 
"open" system. 



C. Asymptotic Variance Formulas 

Here we apply the previous results for statistically homogeneous point patterns to obtain 
asymptotic results for large windows. The conditions under which these expressions yield 
variances that only grow as the surface area of Q are determined. These conditions can be 
expressed in terms of spatial moments of the total correlation function h{r). For simplicity, 
we first consider the case of spherical windows but we show that the results apply as well to 
non-spherical windows. 

Many of our subsequent results will be given for a rf- dimensional spherical window of 
radius R centered at position xq. The window indicator function becomes 

u;(|x-xo|;i?) = e(i?-|x-xo|), (31) 

where Q{x) is the Heaviside step function 

Therefore, the function fi(R), defined in relation (PT|) . becomes the volume of a spherical 
window of radius R given by 

It is convenient to introduce a dimensionless density defined by 

^ = ^''^^'^'>-'wSwf^- '''' 

where D is a characteristic microscopic length scale of the system, e.g., the mean nearest- 
neighbor distance between the points. 

Substitution of the expansion ()A14j) for the scaled intersection volume a{r; R) into ()22|1 
and assuming that the resulting integrals separately converge, yields the variance formula 
for large R as 



(NHR)} - {N(R)y- 



DJ \DJ \DJ 



(35) 
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where A and B are the asymptotic constants given by 

A=l^plhir),r^l^-±j-J^Hr)dr, (36) 

^ = - 2C..(c/2)r(ai)r(i) L ^^'^ 

and i{x) signifies terms of lower order than x j^]. In what follows, the asymptotic constants 
A and B will generically be referred to as "volume" and "surface-area" coefficients for point 
patterns in any dimension. 

Remarks: 

1. Observe that the volume coefficient A is equal to the nonnegative structure factor in the 
limit that the wavenumber approaches zero, i.e., 

A = lim Sik) = l+p f h{r)dr > 0. (38) 

|k|->0 JsRd 

where S(k) is defined by (jl5j) for any dimension. Consistent with our earlier observations 
about relation (|3U|). we see that A is the dominant term for very large windows and indeed 
is the only contribution for infinitely large windows. It is well known that point patterns 
generated from equilibrium molecular systems with a wide class of interaction potentials 
(e.g., hard-sphere, square- well, and Lennard- Jones interactions) yield positive values of A in 
gaseous, liquid, and many solid states. Indeed, A will be positive for any equilibrium system 
possessing a positive compressibility. This class of systems includes correlated equilibrium 
particle systems, an example of which is discussed in Appendix B. The coefficient A will also 
be positive for a wide class of nonequilibrium point patterns. One nonequilibrium example is 



It pa 



To summarize, there is an enormously 



the so-called random sequential addition process 
large class of point patterns in which A is nonzero. 

2. Because the local variance is a strictly positive quantity for i? > [cf. PH]) ]. we have 
from (j35j) that for very large windows 



The crucial point to observe is that if the volume coefficient A identically vanishes, then the 
second term within the brackets of ()35|) dominates, and we have the condition 

B >0, (40) 
14 



where we have used the fact that the variance cannot grow more slowly than ^, i.e., the 
surface area of the window JJj. We will refer to a system in which 

A = lim S(k) = (41) 

as a "hyperuniform" system. Such point patterns do not possess infinite-wavelength fluctu- 
ations. In a recent cosmological study the term "superhomogeneous" has been used to 
describe such systems. Note that for a one-dimensional hyperuniform system, the variance 
is exactly (not asymptotically) given by 

{N\R)) - {N{R)y = 205, (42) 

where B is given by (jTFj) with d = 1, implying that the fluctuations are bounded, i.e., do 
not grow with R j2J]. 

3. By contrast, we will refer to a point pattern in which the surface-area coefficient vanishes 
(-B = 0) as a "hyposurficial" system. A homogeneous Poisson point pattern is a simple 
example of such a system. Inequality ()39|) in conjunction with the fact that the variance 
cannot grow more slowly than the surface area of a spherical (or strictly convex) window 
for statistically homogeneous and isotropic point patterns jlli, enables us to conclude that 
such a system cannot simultaneously be hyperuniform and hyposurficial, i.e., the volume 
coefficient A [cf. (j36j) ] and surface-area coefficient B [cf. (j37|) ] cannot both be zero. In 
Appendix C, we examine the question of how small the volume coefficient A can be made if 
the point pattern is hyposurficial. 

4. Observe also that the asymptotic variance formula (j35j) and the analysis leading to 
condition (pn|) are valid for any statistically homogeneous point pattern. Now if we further 
assume that the point pattern is statistically isotropic, then the volume coefficient ()36p and 
surface-area coefficient (jHTj) can be expressed in terms of certain moments of h, namely, 

A = l + rf2'^0(x'^-i), (43) 

B = - r(^)r(i)^ <- >• 

where 

(x") = / x"/i(x)rfx (45) 
Jo 

is the nth moment of the total correlation function h{x) and x = r/D is a dimensionless 
distance. Following the previous analysis, we see that if A = 0, then the condition for the 
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variance to grow as the surface area implies that the the dth moment of h must be strictly 
negative, i.e., 

(x'^) < 0. (46) 



D. Direct Correlation Function and New Critical Exponents 

The direct correlation function c(r) of a hyperuniform system behaves in an unconven- 
tional manner. In real space, this function is defined by the Ornstein-Zernike equation 

h{r) = c(r) +p f h{r~ r')c{r)dr'. (47) 

This relation has primarily been used to study liquids in equilibrium jl| , but it is a perfectly 
well-defined quantity for general (nonequilibrium) systems, which are of central interest in 
this paper. The second term is a convolution integral and therefore Fourier transforming 
(JUj) leads to 

where c(k) is the Fourier transform of c(r). Using relation (PHj) and definition (jlHjl . we 
can re-express the number variance for a window of arbitrary shape in terms of the Fourier 
transform of the direct correlation function as follows: 



{N'm - {JV(R))2 = {JV(R)) 



1 f5(k;R)_^- 



{27tYJ l-pc(k) 



(49) 



We know that for a hyperuniform system, h{0) = — 1/p by definition, i.e., the volume 
integral of h{r) exists and, in particular, h{r) is a short-ranged function that decays to zero 
faster than |r|~°'. Interestingly, this means that the denominator on the right side of (jlH)) 
vanishes at k = and therefore c(k = 0) diverges to — oo. This implies that the real-space 
direct correlation function c(r) is long-ranged, i.e., decays slower than |r|~°', and hence the 
volume integral of c(r) does not exist. This is unconventional behavior because, in most 
equilibrium instances, c(r) is a short-ranged function, even in the vicinity of thermody- 
namic critical points where h{r) is long-ranged. One can see that c(r) for a hyperuniform 
system behaves similarly to the total correlation function h{r) for an equilibrium system 



near its critical point 



Lriy tc 

la, i 



.e., each of these functions in these respective instances become 
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TABLE I: Definitions of the critical exponents in the vicin- 
ity of or at the hyperuniform state. Here S~^(0) is the 
inverse of the structure factor at A; = 0, ^ is the correlation 
length, and c(r) is the direct correlation function. 
Exponent Asymptotic behavior 

c(r) ~ r'^-d-v (0 = 0^) 



7 

7' 
u 

v' 
V 



long-ranged. If this analogy holds, then one expects the direct correlation function for hype- 
runiform systems to have the following asymptotic behavior for large r = |r| and sufficiently 
large d: 

c(r) ~ {r ^ oo), (50) 

where (2 — rf) < < 2 is a new "critical" exponent associated with c(r) for hyperuniform 



systems that depends on the space dimension 2^. For noninteger values of rj, the asymptotic 
relation (|Hn|l implies that the Fourier transform h(k) is a nonanalytic function of A; = |k|. 
We will show in Section V that there is a class of hyperuniform systems that obey (jKn|) but 
with integer values of t], implying that h(k) is an analytic function of k. Inversion of ()5U|) 
yields 

g(k) ~ {k 0), (51) 

which, when combined with ()48|1 . yields the asymptotic form of the structure factor 

S{k) ~ A;'-" {k 0). (52) 

The specific asymptotic form of 5'(k) for small k contributes to determining the "universal- 
ity" class of the hyperuniform system. 

Let us now consider a point pattern with a reduced density that is nearly hyperuniform 
and which can be made hyperuniform by increasing and/or decreasing the density. We 
denote by 0c the reduced density at the hyperuniform state. The reduced densities and 
0c play the same role as temperature T and critical temperature Tc, respectively, in the 
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analogous thermal problem in the vicinity of a critical point. Thus, we can define critical 
exponents associated with the manner in which certain quantities diverge as the critical 
(hyperuniform) point is approached. For example, for |0c "~ 0| ^ 1, the inverse of the 
structure factor at /c = 0, S'^^(O), and the correlation length ^ obey the power laws 

S-'(O) ~ A-^j ^-C (53) 

f~(l-£)". ^-C. (54) 

where 7 and u are nonnegative critical exponents that are related by the formula 

7 = (2-r/)z/. (55) 

As will be discussed in Section V, ^ characterizes the decay of the direct correlation function 
in the vicinity of = 0c- Analogous critical exponents can be defined for densities near but 
above (pc, as summarized in Table IH In Section VB, we determine the critical exponents 
exactly for certain models of disordered point patterns in d dimensions. 



III. VARIANCE FORMULA FOR A SINGLE POINT PATTERN 

In this section, we derive a new formula for the number variance of a single realization 
of a point pattern consisting of a large number of points in a large system of volume V. 
This is necessarily a volume-average formulation. Fluctuations for a fixed window size arise 
because we let the window uniformly sample the space. As we will show, depending on the 
nature of the point pattern, this formula will generally lead to a result that is different from 
formula (j22I); which was derived for a statistically homogeneous system. We also show that 
the formula derived here is preferable for finding point patterns with an extremal or targeted 
value of the number variance. 

For notational simplicity, we consider a d-dimensional spherical window of radius R, 
keeping in mind that the results of this section apply as well (with obvious notational 
changes) to regular domains of arbitrary shape. We assume that the characteristic size of 
the system is much larger than the window radius so that boundary effects can be neglected 
and that the large numbers ^ 1 and V ^ 1 are comparable such that p = N/V is a finite 
number density. Let us recall relation pn|) for the number of points N{:kq;R) contained 
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within a window at position xq in a system of volume V in which there are points. We 
let the window uniformly sample the space and define the average number of points within 
the window to be 

N 



If" 

^(^) = T7 / E^(|r.-^o|;i?)c?xo 

= p [ e(R-r)dr 
Jv 

= pvi{R) 



D 

where fi(r) and (/) are given by (jHHjl and respectively. 
Similarly, squaring relation (pUj) and averaging yields 

N -I „ N 



(56) 



N^{R) = — / Vw(|ri-xo|;i?)c?xo + — / V iu(|ri - xo|; i?)iy(|rj - xo|; i?)c?xo 



N 



(57) 



where a(r; R) is the scaled intersection volume, given explicitly by ()A5|1 . and 
Therefore, the local variance cr^(i?) is given by 



Ti - r,- 



a'{R) = iV2(i?) - N{R) = N{R) 



N 



l-pv,{R) + -Y.^{T,,-R) 



2V(f)'[l-2^,(^)%lg.(..,fi) 



(58) 



The last term within the brackets is the sum of scaled intersection volumes between all point 
pairs, per point. 

Remarks: 

1. It is important to observe that the series in ()58|l terminates for rij > 2R even for infinitely 
large systems. 

2. Note that the variance formula (jSHjl is different from the ensemble-average formula ()22p . 
which involves an additional weighted average over pairs of points; thus, the appearance of 
the total correlation function h{r). Therefore, the variance function (j58j) . unlike the variance 
function (j22j) . will generally contain small-scale fluctuations with respect to R, of wavelength 
on the order of the mean separation between the points, that are superposed on the large- 
scale variations with respect to R (see examples in Section IV). The expressions and 
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are identically the same for statistically homogeneous (infinite) systems, in which case 
the amplitudes of the small-scale fiuctuations vanish. 

3. Because the variance formula is valid for a single realization, one can use it, in principle, 
to find the particular point pattern which minimizes the variance at a fixed value of R. In 
other words, it is desired to minimize cr'^{R) for a particular value of R among all r^j < 2R, 
i.e., 

min a2(i?) (59) 

Vrij <2R 

where cr'^{R) is given by The scaled intersection volume a(rjj; R) appearing in ()58|) is 

a nonnegative function of Vij (see Fig. Ej) and can be viewed as a repulsive pair potential 
between a point i and a point j. Finding the global minimum of (y^{R) is equivalent to 
determining the ground state for the "potential energy" function represented by the pairwise 
sum in Such global optimization problems can be attacked using simulated annealing 
techniques, for example. More generally, one could devise an optimization scheme in which 
a targeted value of the variance (rather than an extremal value) is sought [26 1. 

4. Because the pairwise sum in fl58|) is positive, we immediately obtain from (j^Hj) the 
following lower bound on the variance: 



a\R) > 2U ( J ) 



This bound is exact for R < where r^^^^^ is the minimum pairwise distance, and 

therefore provides an accurate estimate of the variance for small R. For sufficiently large R, 
however, the bound becomes negative and therefore provides a poor estimate of the variance. 

5. For large R in the special case of hyperuniform systems, the large-scale variations in R 
will grow as R'^~^, and so we have from that 

a\R) = AiR) (g) + O (g) (61) 



where 



d 1 Af 

(62) 



AiR) = 2^ (^) 

is the asymptotic "surface-area" function that contains the small-scale variations in R 



6. It is useful to average the small-scale function A(_R) over R to yield the constant 

A(L) = ^ A{R)dR, (63) 



LJo 
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where A(i?) is given by ()62p. In the case of a statistically homogeneous system, the constant 
surface-area coefficient 

A = lim A(L) = lim ^ A{R)dR (64) 

L-^oo L—*oo L Jo 

is trivially related to the surface-area coefficient B, defined by ()H7|1 in the asymptotic 
ensemble- average formula, by the expression 

A = 2'^4>B = — , y , / h(r)rdr. (65) 



7. Because the formula for the coefficient A is defined for a single realization, we can employ 
it to obtain the particular point pattern that minimizes it. Thus, the optimization problem 
is the following: 

min A, (66) 

all rij<2L 

where A is given by ()63|) . 

8. For large systems in which any point "sees" an environment typical of all points, relation 
()58|) for the variance can be simplified. This requirement is met by all infinite periodic 
lattices for any R as well as statistically homogeneous point patterns for sufficiently large 
R. In such instances, the second term within the brackets of ()58p can be written as sum 
of scaled intersection volumes over — 1 points and some reference point. Thus, we can 
rewrite the variance as 

d N-l 



a\R) 



- 



DJ 



(67) 



fc=i 

where is the distance from the reference point to the fcth point. The asymptotic expression 
(jnU for (T^(_R) and relation for A(_R) still apply but with A(_R) given by the simpler 
formula 



A{R) = 2^0 ( 



Dxd N-1 

We emphasize that the simplified formulas (jU7|) and ()68|) cannot be used for the aforemen- 
tioned optimization calculations. The latter requires the full pairwise sum appearing in the 
general relation (j3Hj) . 

9. In order to make the surface- area function A(/2) or surface-area coefficient A independent 
of the characteristic length scale or, equivalently, density of the hyperuniform point pattern. 
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one can divide each of these quantities by cj)^'^ ^^Z'^, i.e., 



or -JuTTTj,- (69) 



" (j){d-l)/d- 

This scahng arises by recognizing that normahzation of the asymptotic relation (j?)T|l by 
expression for {N{R)) taken to the power (rf — l)/d renders the resuhing normahzed 
relation independent of R/D. Such a scaling will be used to compare calculations of A{R) 
and A for different ordered and disordered point patterns to one another in the subsequent 
sections. Note that since one-dimensional hyperuniform patterns have bounded fluctuations, 
this scaling is irrelevant for d = 1. 



IV. CALCULATIONS FOR INFINITE PERIODIC LATTICES 

It is useful and instructive to compute the variance, using the formulas derived in the 
previous section, for common infinite periodic lattices, which are hyperuniform systems. To 
our knowledge, ex plic it calculations have only been obtained for the square lattice [l^ and 
triangular lattice [1^ in two dimensions. Here we will obtain explicit results for other two- 
dimensional lattices as well as one- and three-dimensional lattices. We take the window to 
be a (i-dimensional sphere of radius R. 

For infinite periodic lattices, Fourier analysis leads to an alternative representation of the 
variance. Let the sites of the lattice be specified by the primitive lattice vector P defined 
by the expression 

P = riiai + n2SL2 H h na-iSLd-i + naSLd, (70) 

where aj are the basis vectors of the unit cell array and rii spans all the integers for i = 
1,2,- ■ - d. Denote by U the unit cell and vc its volume. It is clear that the number of points 
A^(xo; R) within the window at xq [cf. fl20|) ] in this instance becomes 

Ar(xo;i?) = Ee(i?-|P-xo|), (71) 
p 

where the sum is over all P. 

The number A^(xo; R) is a periodic function in the window position xq and therefore it 
can be expanded in a Fourier series as 

iV(xo; R) = pv,{R) + J2 «(q)e''l ' (72) 
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where q is the reciprocal lattice vector such that q ■ P = 27rm (where m = ±1, ±2, ±3 ■ ■ ■) 
and the sum is over all q except q = 0. Following Kendall and Rankin the coefficients 
a(q), for q 7^ 0, are given by 

a(q) = — / A^(xo; i?)e~*^ ' ^orfxo 

Vc Ju 

e(i?-|P-xo|)e-^q-^orfxo 

Vn ^ Ju 



1 

Vc 

1 

Vc J^'^ 

1 / 2vr \ 



Vn J'ik'i 



= '''' 

where Ju{x) is the Bessel function of order v. Note that the integral in the third line is 
nothing more than the Fourier transform of the window indicator function, which is given 
by ()A3|) in Appendix A. The analysis above assumes that there is one point per unit cell, i.e., 
we have considered Bravais lattices. One can easily generalize it to the case of an arbitrary 
number of points nc per unit cell. Formula (fT3|l would then involve nc — 1 additional terms 
of similar form to the original one. 

By Parseval's theorem for Fourier series, the number variance (t'^{R) is given explicitly 

by 

a\R) = — I [N{^o;R)-pvi{R)]^d^o 

Vc Ju 

q^O 

= ^T.(^i]\j^/2mr- (74) 

q^o \ y / 

One can easily obtain an asymptotic expression for the variance for large R by replacing the 
Bessel function in fl74|) by the ffist term of its asymptotic expansion, and thus we have 

/ D \ d—1 / D \ d—2 

2/ T-.\ A / T-.\ / \ . / \ 



a\R) = A{R) (^-j + O (^-j (75) 
where D is a characteristic microscopic length scale, say the lattice spacing, and 



2d+l^d-lj-)2d COS 



qR- {d+l)n/4: 



describes small-scale variations in R. As before, it is convenient to compute the average of 
A(i?) over R to give the surface- area coefficient: 

1 fL 



A = lim - / A{R)dR 

L^oo L Jo 
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^ vl ^^^^ 

It is useful here to apply the specialized volume-average formula (jU7j) to the case of infinite 
periodic lattices. Recognizing that the configuration of an infinite periodic point pattern 
may be characterized by the distances and coordination numbers for the successive 
shells of neighbors {k = 1, 2, 3, . . .) from a lattice point, we find from ()fi7j) that the variance 
can also be represented as 



d r / T3\ d 
1 I _ •}"-rh I 



k=l 



(78) 



The asymptotic expression for <J^{R) and relation (jHSJ for the surface-area coefficient 
A(i?) still apply but with A(i?) given by 



dj I R\ 



A{R) = 2-^0 ( - j 



p s a oo 
k=l 



(79) 



Formula ()78p was obtained by Kendall and Rankin using a more complicated derivation. 
Moreover, their derivation only applies to periodic point patterns. Our more general for- 
mula (jU7j) is also valid for statistically homogeneous point patterns. We also note that our 
most general volume- average representation (j58|) of the variance, from which formula (j67jl is 
derived, is applicable to arbitrary point patterns and its derivation is quite straightforward. 

One can also evaluate the asymptotic coefficient A using the ensemble-average formula 
()65|) . Strictly speaking, this formula is not applicable to periodic point patterns because such 
systems are not statistically homogeneous (neither are they statistically isotropic). To see 
the potential problem that arises by naively applying ()65|) . let the origin be a lattice point in 
the system and consider determining the radial distribution function (72 (r) by counting the 
number of lattice points at a radial distance from the origin. For a lattice in d dimensions, 
we have that 

where Si(r) is the surface area of a sphere of radius r given by p4|) and is the coordination 
number of the kth shell. It is seen that substitution of the corresponding total correlation 
function h{r) = g2{r) — 1 into (jHSI) results in a nonconvergent sum. However, using a 
convergence "trick" 27], one can properly assure a convergent expression by reinterpreting 
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the surface-area coefficient ()65p for a periodic lattice in the following manner: 
A = hm — , ; / e '^'^ h(r)rdr 



lim 



2 



/3^o+ Dr(|; 



e 



(81) 



MDm(3'-^ r(^) 

A. One-Dimensional Examples 

Here we obtain exact expressions for the number of points and number variance for general 
one-dimensional periodic point patterns using the aforementioned Fourier analysis. Using 
this result, we prove that the simple periodic linear array corresponds to the global minimum 
in A. Subsequently, we employ the volume-average and ensemble-average formulations of 
Sections II and III to obtain some of the same results in order to compare the three different 
methods. Recall that hyperuniform systems in one dimension have bounded fluctuations. 



FIG. 3: Portions of two one-dimensional periodic point patterns, where vc = D = 1. The top and 
bottom arrays are the single-scale and two-scale examples, respectively. 

Let us first consider the simplest periodic point pattern in which each point is equi- 
distant from its near neighbors (see Fig. ^ and let this nearest-neighbor distance be unity 
iyc = D = 1). Applying relations (|7^ and (f7H|) and recognizing that q = 2nmai/D 
(m = ±1, ±2, ■ ■ ■) for nonzero q yields that the number of points contained within a one- 
dimensional window of radius R centered at position xq: 

N{R; Xo) = 2R + ^t sin{2nmR) cos{2rrmxo) _ 

According to relation (j74j) . the associated variance is given by 

. 4 I (83) 

The variance cr'^{R) is a periodic function with period 1/2 and is equal to the quadratic 
function 2i?(l — 2R) for < i? < 1/2 (see Fig. H}. Finally, the surface- area coefficient A, 

25 



0.4 



0.3 



DC 



0.2 



0.1 



1 1 1 1 1 1 1 

Single-Scale Linear Lattice 









DC 



0.5 



1 

R 



1.5 




FIG. 4: Left panel: The quadratic periodic variance function cr'^{R) for the single-scale periodic 
one-dimensional point pattern given by (jHHI)- The horizontal line is the average A = 1/6. Right 
panel: The piecewise-quadratic periodic variance function o"^(i?) for the two-scale periodic one- 
dimensional point pattern given by ()87() for the case C = 1/4. The horizontal line is the average 
A = 7/24. 

defined by (|77p. which in one dimension amounts to the positional average of the variance 
for any value of R, is exactly given by the constant 



1 OO 1 



m=l 



m 



_ 1 
2 - g- 



(84) 



It is known that this simple linear array yields the minimum value of A among all one- 
dimensional regular lattices. This is intuitively clear from the volume-average variance 
relation (j^Hj) for d = 1; the linear repulsive effective "pair potential" contained therein is 
evidently responsible for such a minimum. However, heretofore it was not known whether 
this pattern corresponded to a global minimum, i.e., the smallest value of A among all 
infinite one- dimensional hyperuniform patterns. We now prove that the single-scale lattice 
indeed produces the global minimum. To prove this assertion, we utilize the identity 



TT 



1^1 + 2 cos(27rmx) 1 
1 L — = o - 2a;(l - x) 



m=l 



(85) 



and note that f{x) is a convex quadratic nonnegative function for all real x. Now consider 
a case in which there are M points per unit cell in which the length of the unit cell is still 
unity. Thus, excluding the point at each lattice site, there are M — 1 points inside the unit 
cell with positions Ci, C2, • • • , Cm-i such that each C,i lies in the interval (0, 1). Without loss 
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of generality, we arrange the M — 1 points such that Q < Q+i {i = 1, 2, . . . , M — 2), but their 
positions are otherwise arbitrary. Following a similar analysis as the one above, we find that 
the number of points within a window centered at Xq is exactly given by 

°° sin(27rmi?) f Ef-n^ cos[27rm(xo - d)] 
N{R; xo) = 2MR + 2 ^ — '-^^ l^JllHi (gg) 

m=l 

where Co = 0- The variance is therefore given by 

2 ^ sm\27TmR) \m + cos(27rmO) + E*<fc' cos[27rm(a - 0)]! 

^ = — E 2 -■ 87 

We see that the variance (t'^{R) for an arbitrary one-dimensional point pattern within the 
unit cell is a periodic function with period 1/2. (As we will see, the variance in higher 
dimensions is not a periodic function in R for periodic point patterns.) The average of the 
variance is exactly equal to the surface-area coefficient ()77|) : 

_ ^ ^ - M + Ejll' cos(27rmO) + cos[27rm(Cfc - Q)] 

(88) 

j=i j<k 

where f{x) is given by ()85|) . Because A is given by a sum of convex quadratic nonnegative 
functions, the global minimum is found from the zeroes of the derivative dA/d(n'- 

dA ""^ ^^^^ 

— = = l-2C„ + E(l-2Cn + 20)- E (l + 2Cn-20), (ri = l,2,...,M-l). (89) 

'^Sn. j=i j=n+l 

It is easy to verify that the global minimum is achieved when the M — 1 are uniformly 
distributed in the interval (0, 1), i.e., (n = n/M {n = 1, 2, . . . , M — 1), yielding A = 1/6. 
Since this result is valid for arbitrary M, the simple single-scale lattice produces the global 
minimum value of A among all infinite one-dimensional hyperuniform point patterns. 

Note that the single-scale lattice corresponds to the densest packing of one-dimensional 
congruent hard spheres (rods) on the real line. This might lead one to conjecture that the 
Bravais lattice associated with the densest packing of congruent spheres in any space dimen- 
sion d provides the minimal value of A among all periodic lattices for spherical windows. As 
we will see, this turns out to be the case for d = 2, but not for d = 3. 

The variance as computed from (jH7|) for the case M = 2, which we call the "two-scale" 
lattice (see Fig. E)), is included in Fig. IHfor C = Ci = 1/4. In this instance, A = 7/24. 
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TABLE II: The surface- area coefficient A for some ordered and disordered one-dimensional point 
patterns. The result for the two-scale lattice is for = = 0.25. 



Pattern 


T 


A 


Single-Scale Lattice 


1 


1/6 K, 0.166667 


Step-|-Delta-Function §2 


0.75 


3/16 = 0.1875 


Step-Function §2 


0.5 


1/4 = 0.25 


Two-Scale Lattice 


2 


7/24 K, 0.291667 


Lattice-Gas 


1 


1/3 f« 0.333333 



Clearly, the variance for the two-scale case bounds from above the variance for the single- 
scale case. Table 101 compares the surface-area coefficient for the single-scale and two-scale 
one-dimensional lattices. The other one-dimensional results summarized in Table ITTl will be 
discussed in the ensuing sections. The potential use of the local variance as an order metric 
for hyperuniform point patterns in any dimension is discussed in Section VI. 

Consider obtaining the volume-average representation of the variance for the two afore- 
mentioned one- dimensional periodic patterns from ()78j) . Using relations ()A5|1 and ()A9|1 
for the scaled intersection volume a(r;i?), we find for any one-dimensional periodic point 
pattern in which D = 1 that 

Mr 



a\R) = 2(f)R 



1 - 20i? + ^ Zfc ( 1 - H j 0(2i? - 



(90) 



k=l 

where Mfj corresponds to the largest value of k for which < 2R. Because in one dimension 
A(_R) = ct'^{R), where A(_R) is the function defined by ^B^, it follows that the average A is 
given by 



A^2f\(i.).i.^,4(l-^)-,E^.(l-^/ 



(91) 



40L _^ 

' k=l 

where Ml corresponds to the largest value of k for which < 2L. Using the fact that 
(j) = 1, Tk = k, and Z/^ = 2 for all k for the single-scale lattice, one can easily reproduce the 
graph for cr^(i?) depicted in Fig. HI using relation and verify that A = 1/6 employing 
relation (jH^)- Similarly, for the two-scale case, we have that (p = 2, Vk = k/A and Zk = 1 for 
odd k, and = k/2 and Z^ = 2 for even k. Hence, relation ()90|) leads to the same graph of 
the variance shown in Fig. EJ and relation yields A = 7/24 for ( = 0.25, as before. 
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We can also compute the surface-area coefficient using the ensemble-average relation ()8ip. 
In one dimension, this relation yields 



A 



lim 



k=l 



(92) 



where we have taken D = 1. The sum in can be evaluated exactly using the Euler- 
Maclaurin summation formula |23]. If f{k) is a function defined on the integers, and con- 
tinuous and differentiable in between, the Euler-Maclaurin summation formula provides an 
asymptotic expansion of the sum J2k=ofi^) as ^ oo. Applying this asymptotic formula 
to (jn^ in the cases of the single-scale and two-scale lattices, yields that A = 1/6 and 
A = 7/24, respectively, which agree with the results obtained using the previous two meth- 
ods. Although the Fourier-analysis and volume-average procedures are more direct methods 
to determine A for one- dimensional lattices, we will see that the representation ()8H) provides 
an efficient means of computing A for lattices in higher dimensions. 

1 
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FIG. 5: The asymptotic surface-area function A(i?) for the square lattice for 1 < < 4, where D 
is the lattice spacing. The horizontal line is the asymptotic average value A = 0.457649. 



B. Two-Dimensional Examples 

Here we evaluate variance characteristics for the following four common two-dimensional 
lattices: square, triangular, honeycomb, and Kagome lattices. From the lattice series (ffljl . 
fl7fj|l and ()77j) with d = 2, we have general two-dimensional series relations for the variance 
a^{R), asymptotic surface-area function A(i?) and surface-area coefficient A, respectively. 
For a specific lattice, the evaluation of any of these series requires the reciprocal lattice 
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vector q and Vc- For example, for the square lattice, q = 27r[miai + m2SL2\/D {rrii = 
0, ±1, ±2, • ■ ■) for nonzero q and vc = . The sums are straightforward to evaluate, even if 
they converge slowly. Provided that R is not very large, however, the corresponding volume- 
average relations (j78|) and (j79|) are superior for computational purposes because the series 
involved are finite rather than infinite. For example, the asymptotic surface-area function 
A(i?) for the square lattice is plotted in Fig. El using (fTUI) with ()A10|) for 1 < i? < 4. The 
function is seen to be aperiodic, but fluctuates around an average value in a bounded fashion. 
It is noteworthy that the behavior of A(i?) for larger values of R is qualitatively the same. 
Interestingly, the average value of A(i?) over this small interval near i? = (as well as other 
intervals of the same length) is quite close to the infinite- interval average value A [2^. 

The average value of the surface-area function A(i?) over all i?, equal to the surface- area 
coefficient A [cf. (fTTj) ]. is given (to six significant figures) by A = 0.457649. The series (f77|) 
for the square lattice was first evaluated by Kendall Because it is a slowly converging 
series, he exploited certain results of number theory to re-express the sum in terms of a more 
rapidly convergent series. 

TABLE III: The surface-area coefficient A for some ordered and disordered two-dimensional point 
patterns. For ordered lattices, (j) represents the close-packed covering fraction. 



Pattern 




A/(Ai/2 


Triangular Lattice 


7r/^/l2 « 0.9069 


0.508347 


Square Lattice 


7r/4 « 0.7854 


0.516401 


Honeycomb Lattice 


7r/(3^/3) « 0.6046 


0.567026 


Kagome Lattice 


37r/(8^/3) K, 0.6802 


0.586990 


Step-|-Delta-Punction 52 


0.5 


25/V(3vr) w 0.600211 


Step-Function g2 


0.25 


8/(37r) K, 0.848826 


One-Component Plasma 




1.12838 



We found that numerical evaluation of the ensemble- average relation ()8H1 is a simple 
and effective means of computing accurately the surface-area coefficient A for any common 
lattice. In two dimensions, this relation yields 



A 



lim 



160^ 80 - ^ - 
^1/2^3/2 ^ .tl "" 



(93) 
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The sum in fl93|) is easily computed as a function of the convergence parameter (3 for any 
simple lattice. For sufficiently small (3, this sum is linear in (3 and extrapolation to /3 — ^ O"*" 
yields results that are accurate to at least six significant figures. We have also computed 
the surface-area coefficient for triangular, honeycomb, and Kagome lattices. The result for 
the triangular lattice was first reported by Kendall and Rankin jl^, but the results for 
the honeycomb and Kagome lattices are new. In Table IIIH we compare all of these re- 
sults for the common two-dimensional lattices to one another by tabulating the normalized 
.a...depende. s.face-a.ea coefflcie,., ,.e., A/,- |cf. 63], RanMn proved that 
the triangular lattice has the smallest normalized surface-area coefficient for circular win- 
dows among all infinite periodic two-dimensional lattices, which is borne out in Table 11111 
However, there is no proof that the triangular lattice minimizes among all infinite 

two-dimensional hyperuniform point patterns for circular windows. Included in Table 
are results for disordered point patterns that will be discussed in the ensuing sections. 

1 



0.5 



o 

J" 
\_ 



-0.5 





100 



102 



104 106 
R/D 



108 



110 



FIG. 6: The difference between the normalized scale-independent surface-area function A(i?) for 
the triangular and square lattices as a function of R, where D is the lattice spacing. Here A^j.j = 
A(i?)tri/0;/- and Asq = A(i?)sq/,/.^^l 



Although the normalized surface- area coefficient is smallest for the triangular lattice, 
Table ITTTl reveals that the corresponding coefficients for the other lattices are not appreciably 
larger. This suggests that the fiuctuating surface-area function A(i?) for non-triangular 
lattices may be smaller than the corresponding function for the triangular lattice for certain 
values of R. This is indeed the case as illustrated in Figure IHl where the difference between 
the normalized scale-independent surface-area function A(i?) for the triangular and square 
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lattices is plotted for the range lOOD < R < llOD using relation ()79|) . This difference 
oscillates rapidly about zero over this range of R, but the same qualitative trends occur for 
all values of R and for any pair of periodic lattices considered here. Given our previous 
interpretation of the global minimum of the variance as corresponding to the ground state 
of a many-particle system with a potential energy function given by a{r; R) (Section III), we 
see that the optimal lattice structure is sensitive to small changes in the value of R (which 
determines the range of the potential). This calls into question previous studies 3l| that 
claim to have found stable ground-state lattices for two-dimensional systems of particles with 
purely repulsive interaction potentials of the same qualitative form as shown for a{r; R) in 
Fig. Hwith d = 2. 



C. Three-Dimensional Examples 

Here we specialize to common infinite three-dimensional periodic lattices: simple cubic 
(SC) lattice, face-centered cubic (FCC) lattice, hexagonal-close-packed (HCP) lattice, body- 
centered cubic (BCC), and the diamond lattice. Explicit results for the number variance 
for such lattices have heretofore not been reported. From (jZH), ()7(i|l and (fTTjl with d = 3, 
we have general three-dimensional relations for the variance a'^{R), asymptotic surface- area 
function A{R) and surface-area coefficient A, respectively. These expressions are easily 
evaluated for the specific lattice given vc and the reciprocal lattice vectors q. As we noted 
earlier, the volume-average relations ()78|) and ()79p for c/ = 3 are superior for computational 
purposes provided that R is not very large. Qualitatively, the three-dimensional trends for 
the surface-area function A{R) are similar to the two-dimensional ones described above (see, 
for example. Figs. El and ^ and so we will not explicitly present such three-dimensional 
results here. 

The ensemble-average relation ()81|) . which for d = 3 and D = 1 yields 



A = lim 

(3^0+ 



(94) 



and provides an efficient means of computing the surface-area coefficient A for three- 
dimensional infinite periodic lattices by extrapolating the results for sufficiently small f3 to 
(3^0^ This has been carried out for all of the aforementioned common three-dimensional 
lattices and the results are summarized in Table IIV[ where we tabulate the normalized 
scale-independent surface-area coefficient, i.e., A/^^/s [cf. fp|) ]. 
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TABLE IV: The surface-area coefficient A for some ordered and disordered three-dimensional 
point patterns. For ordered lattices, (p represents the close-packed covering fraction. 



JT cxL IjCl 11 






BCC Lattice 


07r/(8v,ij U.bsUz 


1.24476 


1^ Lattice 


vr/ V i» ~ U.74UO 


1.24552 




7r/\/T8 « n 7AV\'\ 




SC Lattice 


7r/6 « 0.5236 


1.28920 


Diamond Lattice 


37r/(16^/3) « 0.3801 


1.41892 


Damped-Oscillating gi 


0.46 


1.44837 


Step+Delta-Function g2 


0.3125 


5I/3 . 9/210/3 p. 1.52686 


Step-Function g2 


0.125 


2.25 



Contrary to the expectation A/(f)'^^^'^~^^ should, among all lattices, be a global minimum 
for the closest-packed lattices for spherical windows, we find that the minimum in three 
dimensions is achieved for the BCC lattice, albeit very close in numerical value to the FCC 
value (the next smallest value) j^^l- This suggests that the closest-packed Bravais lattice 
for (i > 3 does not minimize A/(f)'^^^'^~^^ 2^. Included in Table HVl are results for disordered 
point patterns that will be discussed in the ensuing sections. 



V. NON-PERIODIC HYPERUNIFORM SYSTEMS 



In this section, we briefiy describe the known non-periodic hyperuniform point patterns 
in one, two, and three dimensions and identify some others. For certain one-, two, and three- 
dimensional disordered hyperuniform point patterns, we exactly determine the corresponding 
surface-area coefficients, structure factors, direct correlation functions, and their associated 
critical exponents. A discussion concerning the potential use of surface-area coefficient A as 
an order metric for general hyperuniform point patterns is reserved for Section VI. 
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A. Examples 



Statistically homogeneous hyperuniform point patterns in one dimension are not difficult 
to construct. Two examples are discussed here: one is a "lattice-gas" type model and the 
other is a construction due to Goldstein et al. js^. The first example is constructed by 
tesselating the real line into regular intervals of length D. Then a single point is placed 
in each interval (independently of the others) at any real position with uniform random 
distribution. The number density p = 1/D, and the pair correlation function is simply given 
by 

f r/D, r < D, 
92{r) = { ' ' - ' (95) 

One can easily verify that the system is hyperuniform [A = 0) and that the surface-area 
coefficient (jU^ is given by 

A = ^, (96) 

exactly twice the surface-area coefficient for the simple single-scale periodic point pattern [cf. 
(jHll)]. This one-dimensional lattice-gas model is a special case of the so-called d-dimensional 
shuffied lattice that we will describe below. 

A less trivial example of a statistically homogeneous one-dimensional hyperuniform sys- 
tem is the construction of Goldstein et al. 3^, which obtains from a homogeneous Pois- 
son point process a new hyperuniform point process. This construction is defined as fol- 
lows: First, one defines a statistically homogeneous process X{x) on the real line such that 
X{x) < 1. This process is specified by dynamics such that X{x) decreases at the rate of 
unity, except at the points of the Poisson process, where X{x) jumps up by one unit unless 
this jump violates the upper bound condition, in which case no jump occurs. Second, one 
takes the points of the new point process to be those points in which X{x) actually jumps. 
This new point process is hyperuniform. It is not known how to extend this construction to 
higher dimensions {d > 2). 

The construction of statistically homogeneous and isotropic point patterns that are hy- 
peruniform in two or higher dimensions is a challenging task. An example of a statisti cally 
homogeneous (i-dimensional system that is hyperuniform is the so-called shuffled lattice j35| , 
but it is not statistically isotropic. This is a lattice whose sites are independently randomly 
displaced by a distance x in all directions according to some distribution with a finite second 
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moment. 

Gabrielli et al. 35 1 have observed that the a point pattern derived from the "pinwheel" 
tihng of the plane 36| has a number variance that grows as the surface-area (perimeter) of 
the window, and is statistically homogeneous and isotropic. The prototile of the pinwheel 
tiling is a right triangle with sides of length one, two, and \/5. The tiling is generated by 
performing certain "decomposition" and "inflation" operations on the prototile. In the first 
step, the prototile is subdivided into five copies of itself and then these new triangles are 
expanded to the size of the original triangle. These decomposition and inflation operations 
are repeated ad infinitum until the triangles completely cover the plane (see Fig. Ej). It is 
obvious from the aforementioned discussion that the point pattern that results by randomly 
placing a point in each elementary triangle is hyperuniform. Importantly, because the tiles 
appear in infinitely many orientations, one can show that the resulting pattern is not only 
statistically homogeneous but statistically isotropic. The full rotational invariance of the 
pattern is experimentally manifested by a diffraction pattern consisting of uniform rings 
rather than isolated Bragg peaks. 




FIG. 7: Portion of a pinwheel tiling. 



The one-component plasma is a statistical mechanical model that is known to have a num- 
ber variance that grows only as the surface area of the window [17|, |37[ . The one-component 
plasma is a system of point particles of charge e embedded in a uniform background that im- 
parts overall charge neutrality. In d = 2, the n-particle correlation functions for this model 
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are exactly solvable in the thermodynamic limit when the coupling constant F = e^/ (kT) = 2 
38[, and, in particular, the total correlation function is then given by 



h{r) = -e""'"' . (97) 



Substitution of into (jHTj) gives the surface area coefficient [37| as 



A = ^0^/^ (98) 



where = pnD^/A. This evaluation of A is included in Table UTTl Observe that the structure 
factor of the (i-dimensional one-component plasma at small k behaves as 

S{k) {k^ 0) (99) 

and, therefore, the corresponding asymptotic behavior of the Fourier transform of the direct 
correlation function is given by 

c(A;)~-^ (A;^0). (100) 



Another interesting model that is known to be hyperuniform 



Zeldovich 



35l | is the Harrison- 



39[ power spectrum for the primordial density fluctuations in the universe. Here 



the structure factor for small k behaves as 

S{k) ~ k. (101) 



Recently, Gabrielli et al. [35] have discussed the construction of point patterns in three 
dimensions that are consistent with the Harrison- Zeldovich spectrum. 

The pre sent authors have recently introduced and studied so-called (72-iiivariant processes 
40,1^|. A (72-invariant process is one in which a chosen nonnegative form for the pair 
correlation function g2 remains invariant over a nonvanishing density range while keeping 
all other relevant macroscopic variables fixed. The upper limiting "terminal" density is the 
point above which the nonnegativity condition on the structure factor [cf. (jl8|) ] would be 
violated. Thus, at the terminal or critical density, the system is hyperuniform if realizable. 
In the subsequent subsection, we will calculate the surface-area coefficient exactly for several 
of these (yf2-invariant processes. We will also exactly determine the corresponding structure 
factors, direct correlation functions, and their associated critical exponents. 
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Interes ting ly, random packings of spheres near the maximally random jammed (MR J) 



state j4^, |43| appear to be hyperuniform. Figure |H1 depicts the structure factor for such 
a computer-generated 40,000-particle packing, is vanishingly small for small wavenumbers. 
The packing is strictly jammed j^l, which means that the particle system remains me- 
chanically rigid under attempted global deformations (including shear) that do not increase 
volume and, furthermore the packing is saturated. A saturated packing of hard spheres 
is one in which there is no space available to add another sphere. In the case of satu- 
rated packings of identical hard spheres of unit diameter, no point in space has distance 
greater than unity from the center of some sphere. An attractive postulate would be that 
all strictly jammed saturated infinite packings of identical spheres are hyperuniform. Ex- 
amples of strictly jammed saturated periodic packings in two and three dimensions include 
the closest packed triangular and face-centered cubic lattices, respectively. In light of this 
discussion, one can view a disordered packing near the MRJ state as a type of "glass" for the 
hard-sphere system. An important open fundamental question is whether there are molec- 
ular glasses (with "soft" intermolecular potentials) that become hyperuniform in the limit 
that the temperature vanishes. Indeed, our preliminary results indicate that this possibility 
is attainable. 
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FIG. 8: The structure factor for a random packing of three-dimensional identical hard spheres of 



diameter D near the MRJ state 
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as computed from a single realization consisting of 40,000 
^articles in a cubical box with periodic boundary conditions using the protocol described in Ref. 
4^ . The packing (covering) fraction of spheres (/> is 0.632. 
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B. Exact Results for (72-Invariant Processes 



Here we evaluate the surface-area coefficient exactly for three different disordered g2 



invariant processes studied by us earlier 



20 



4l|. We also exactly determine the corre- 



sponding structure factors, direct correlation functions, and their associated critical expo- 
nents. 



1. Step-Function §2 

Let us first consider the (yf2-invariant process in which a spherically symmetric pair corre- 
lation or radial distribution function is defined by the unit step function [4^ 

f 0, r<D, 

92{r) = Q{r-D) = { " (102) 

[1, r>D, 

The condition g2{r) = for r < D prevents any pair of points from getting closer than 
a distance D to one another. Note that in the special case of a system of identical hard 
spheres in equilibrium in the limit p ^ 0, g2 is exactly given by (|102j) . The corresponding 
total correlation function is given by 

f -1, r<D, 
h{r) = -e{D - r) = { - (103) 

[0, r>D, 

which when substituted into and yields the volume and surface-area coefficients as 
A = Sik = 0) = l- 2\ B = ^ = -JL:^fI^^4^$—(f). (104) 

^ ^ ^' 2<i^ r((rf + 3)/2)r(i/2)^ ^ ^ 

The reduced density defined by (equivalent to the covering fraction of the hard cores 
of diameter D) lies in the range < (p < (pc, where 

^c=Y. (105) 

is the terminal or critical density, i.e., the density at which the system is hyperuniform, 
where ^ = and 

d'Tid/2) 

^ = ^=4r((rf + 3)/2)r(i/2)- 

The values of the scale-independent surface-area coefficient A/0('^-^)/'^ for d = 1,2 and 3 
are given in Tables II, III and IV, respectively. It is noteworthy that a recent study 
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provides convincing numerical evidence that the step-function g2 is reahzable by systems of 
impenetrable dimensional spheres (with d = 1 and d = 2) for densities up to the terminal 
density. Thus, it appears that satisfying the nonnegativity conditions on g2ir) and S{k) in 
this instance is sufficient to ensure realizability. 

The Fourier transform of the total correlation function (|1(J3|) yields the analytic function 

/ 27r \ '^/'^ 

h{k) = - (— j D'j,„{kD). (107) 
Thus, use of flT3j) gives the structure factor for (p in the range < < 0c to be 

s{k) = 1 - r(i + d/2) (^)'^' [£j JctnikD). (108) 

Similarly, the Ornstein-Zernike relation ()48|1 yields an exact expression for the Fourier trans- 
form of the direct correlation function: 



/ 27r \ "'/^ 



m = ^ ■ (109) 



Thus, the small-A; expansions of S{k) and c{k), which determine their behavior in the 
vicinity of, and at, the critical point, are respectively given by 

^«=(^-|) + 2(^£<*^)' + '^'<*^''' 

and 

where fi(-D) is the volume of a (i-dimensional sphere of radius D [cf. fl33|) ]. At the critical 
point = 0c, we see that S{k) ~ k'^ and c{k) ~ —k'"^, and therefore comparison to (|^ 
and yields the exponent r] = 0. Relation (jllO|) leads to the power law 

S-\0)=(^1-^^ , 0^0c-, (112) 

which upon comparison to (j^Hj) immediately yields the critical exponent 7=1. The corre- 
lation length ^ is defined via (|llip . which we rewrite as 

k^c{k) + r^d{k) = -G, kD < 1 (113) 
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where 



D 



-1/2 



G 



[2(rf + 2)0jV2 
2{d + 2)vi{D) 



1 - — 



(114) 

(115) 

to the 



and vi{D) is the volume of a sphere of radius D defined by (jH^ . Comparison of 
power law yields the exponent u = 1/2. Note that the exponent values 7 = 1,^ = 1/2, 
and 1] = are consistent with the interrelation (j55j) . Inversion of (jll3|) yields the partial 
differential equation 



V^c(r) - r^c(r) = GS{r), r :^ D, 



(116) 



where the spherically symmetric Laplacian operator in any dimension d is given by 



1 d 



dr 



We see that the direct correlation function in real space for large r is determined by the 
Green's function of the linearized Poisson-Boltzmann equation. 

Let us first determine the solutions of ()116j) at the critical point (j) = (f)c where ^ diverges 
to infinity. Thus, the asymptotic behavior of c(r) for r ^ D is given by the infinite-space 
Green's function for the (i-dimensional Laplace equation j42], and so we obtain 

d=l, 

d = 2, 



c r 



Ak 

Hi. 

2{d+2) fr^ ^^-2 



(118) 



d> 3. 



d{d-2) 

Observe that it is only for d > 3 that c(r) follows the power-law form (j51|) with an exponent 
7] = 0. The fact that 77 takes on an integer value is due to the fact that h{k) is an analytic 
function. Note also that the real-space direct correlation function of the one-component 
plasma has precisely the same asymptotic form as pi8p . albeit with different amplitudes 
(prefactors). 



As cxD for fixed r, the solutions of pi6|l are 



c(r) 



(I) ^^P(-^/^)' 



4— In ( — ) exp(-r/0. 



d=l, 
d = 2. 



(119) 



2{d + 2)<P, fry-' I ... 
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On the other hand, it is noteworthy that as r — oo for fixed ^, the asymptotic behavior 
changes according to the relation 

2. Step + Delta Function §2 

Here we consider the 5f2-invariant process defined by a radial distribution function that 
consists of the aforementioned unit step function plus a delta function contribution that acts 
aX r = D: 

g^[r) = e(r -D) + -^-^5ir - D), (121) 
pSi{D) 

where Z is a nonnegative constant and Si{D) is the surface area of a sphere of radius 
D defined by 
and Stillinger 



f|). The function p2H) was one of several examples studied by Torquato 
I to understand the relationship between short-range order and maximal 
density in sphere packings. In this investigation, Z was interpreted as the average contact 
coordination number. Here we consider their case IV (given in the appendix of Ref. j2o|) in 
which the condition 

Z = 4^<P (122) 



is obeyed in order to constrain the location of the minimum of the structure factor to be at 
/c = 0. Here the reduced density lies in the range < < 0c, and 

rf + 2 , ^ 

^^ = ^ (123) 

is the terminal or critical density. Note that the function specified by relation ^V2\\ji is a 
special limit of the radial distribution function corresponding to the dilute and narrow limit 
of the square-well potential studied by Sakai, Stillinger, and Torquato 

Substitution of ()121|) into (jlSj) and (jl^ yields the volume and surface-area coefficients as 

A = S(k = 0) = ,-—4>. ^-^= (,^2)r((d + 3)/2)r(l/2) ^- 
At the critical density, A = and 

, d\d+2nd/2) 

^ = ^^^B= 16r((rf + 3)/2)r(l/2)- ^^25) 
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The values of the scale-independent surface-area coefficient A/(j)^'^~^^/'^ for d = 1,2 and 3 are 
given in Tables HH IIIII and IIV| respectively. 

The combination of relations ()15|) . ()48|) . and p21|) give the structure factor and Fourier 
transform of the direct correlation function respectively for (p in the range < (p < (pc- 



S{k) 



2'^/'^T{2 + d/2) / 

+ (A;D)W2)-1 [j^ 



J(d/2)-i{kD) Jd/2{kD) 



d + 2 



kD 



{27iY/^D^ 


J{d/2)-i{kD) Jd/2{kD) 
d + 2 kD 




2'^/^r{2 + d/2) 
^ (fcD)W2)-i 


\<Pc) 




' J^a/2)-i{kD) Jd/2{kD) 
d + 2 kD 



(126) 



(127) 



Therefore, the Taylor expansions of S{k) and c{k) about A; = are respectively given by 

1 



and 



S{k) 



1-- + 



i{d + 2){d + 4) 0, 
-2v,{D) 



ikDY + 0[{kDf 



128) 



(129) 



Relation 



V (pc) 8(rf + 2)(rf + 4)0, 
leads to the power law 

/ , \ -1 



ikDf + 0[{kDf 



s-\o) 



1 - — 



(130) 



which upon comparison to ()53|) again yields the critical exponent 7 = 1. The correlation 
length ^ is defined via p29|) . which we rewrite as 



where 



G 



k%k) + r c{k) = 

D 

[8(ci + 2)(ci + 4)0,]i/4 
16{d + 2){d + 4:)vi{D) 



kD < 1 



-1/4 



1 - 



(131) 

(132) 
(133) 



D^ 

Comparison of p32j) to the power law (jH^ yields the exponent u = 1/4. We see that the 
exponent values 7 = 1, ^ = 1/4, and 77 = —2 are consistent with the interrelation ()55j] . 
Inversion of ()131|) yields the partial differential equation 



V^c(r) + r^c(r) = -G6{r), r :^ D, 



(134) 
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where = V^V^ is the spherically symmetric biharmonic operator, and is given by 

(cnD. 

The solutions of (jl34j) at the critical point = 0c — ^ oo) are given by the infinite-space 
Green's function for the rf-dimensional biharmonic equation. It is only for d > 5 that the 
solutions admit a power law of the form (j54j) with an exponent f] = —2, namely, 

3. Damped- Oscillating g2 



In three dimensions, Torquato and Stillinger j20| also considered a (yf2-invariant process 
that appends a damped-oscillating contribution to the aforementioned step+delta function 
g2- Specifically, they examined the radial distribution function 



p4"^ 



g^{r) = e(r -D)+ (5(r - D) + — e-"^'' sin(a3r + a4)e(r - D). (136) 



Here we consider their case II, where at the terminal density 0c = 0.46, Z = 2.3964, 
Oi = 1.15, 02 = 0.510, 03 = 5.90, and 04 = 1.66. At this critical point, the volume coefficient 
A = and the surface-area coefficient is given by 



A = 360^ - 60cZ + 144ai0^/, (137) 



where 



/oo 

{2al — a| — 60303 — 20303 — 402O3 — 402O3 — O2O3) 



e "''^ cos(o3 + 04) 



3; 



(2O3 + 6O2O3 — O2 — 2O2 — 2O2 — O2O3 — 2020^ 

(oi + oi 



f^-^ — ,,2 , 2N3 — ^e-"2 sin(o3 + 04). 



Substitution of the aforementioned parameters in ()137p yields A = 0.863082. This evaluation 
of A is included in Table HVl With this choice of (72, the first non-zero term of the small-fc 
expansion of the structure factor S{k) at the critical point is of order /c"^, and therefore the 
exponent r] = —2, as in the previous case. However, here c(r) does not admit the power-law 
form (50) for large r because rj < —1. 
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VI. DISCUSSION AND CONCLUSIONS 



The principal theme presented in this paper is that number fluctuations calculated for 
variable window geometries offer a powerful tool to characterize and to classify point-particle 
media. This theme encompasses both spatially periodic (crystalline) particle patterns, as 
well as those that are globally disordered (amorphous). By considering the large- window 
asymptotic limit, special attention attaches to "volume" and to "surface" fluctuations in 
space dimension d>l. A special class of "hyper uniform" point patterns has been recognized 
for which the volume fluctuations vanish identically; equivalently these are systems for which 
the structure factor S{k) vanishes at A; = 0. Another special class of "hyposurficial" point 
patterns has also been recognized for which the surface fluctuations vanish identically. The 
first of these special attributes requires that the {d — l)-st spatial moment of the total 
correlation function be constrained in magnitude; the second requires a similar constraint on 
the d-th spatial moment of the total correlation function. The preceding text demonstrates 
that no point pattern can simultaneously be both hyperuniform and hyposurficial. 

All infinitely extended perfectly periodic structures are hyperuniform. We have stressed 
that geometrically less regular cases of hyperuniformity also exist, including those that 
are spatially uniform and isotropic. The suitably normalized surface fiuctuation quantity, 
which measures the extent to which hyperuniform systems fail to attain hyposurficial status, 
becomes a natural nonnegative order metric that we have evaluated numerically for a basic 
sampling of structures. We proved that the simple periodic linear array yields the global 
minimum value for hyperuniform patterns m d = 1, and showed that the triangular lattice 
produces the smallest values for the cases tested in ci = 2. But in spite of the fact that 
these minimizing structures correspond to optimal packings of rods and disks, respectively, 
the face-centered-cubic lattice for optimal sphere packing does not minimize the surface- 
fiuctuation order metric for ci = 3. Instead, the body-centered cubic lattice enjoys this 
distinction For each choice of space dimension, other lattices and irregular hyperuniform 
patterns yield higher values for this order metric. An order metric for hyperuniform systems 
based on the local variance may find potential use in categorizing "jammed" and "saturated" 



sphere packings 42, |43|, |4^ |47[ whose long- wavelength density fiuctuations vanish. 



It is clearly desirable to extend the set of point patterns for which the surface fiuctuation 
order metric has been numerically evaluated. This would help to strengthen the impression 
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created thus far that regardless of space dimension d, point patterns arranged by increasing 
values of the order metric are indeed essentially arranged by increasing structural disorder. 
It will be important in the future to include a selection of two and three-dimensional qua- 
sicrystalline point patterns 48] in the comparisons; the presumption at the present state of 
understanding is that they would present order metrics with values that lie between the low 
magnitudes of periodic lattices, and the substantially larger magnitudes of spatially uniform, 
isotropic, irregular point patterns. It would also benefit insight to include cases of spatially 
uniform, but anisotropic, point patterns; for example, those associated with "hexatic" order 
in two dimensions |49|. 



An irnportant c 
cesses" |20 



4G 
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ass of hyperuniform systems arises from the so-called "(yf2-invariant pro- 



45|. These processes require that the pair correlation function g2{r) 



remain unchanged as density increases from zero. For those 5f2-invariant processes that cor- 
respond to thermal equilibrium, this criterion is implemented by virtue of compensating 
continuous changes in the particle pair potential function. For any given choice of the in- 
variant g2, such a process is in fact achievable, but only for densities up to a terminal density 
limit. At this upper limit, the system of points attains hyperuniformity, i.e., S{k) = 0. Fur- 
thermore, examination of the Ornstein-Zernike relation reveals that the direct correlation 
function c(r) develops a long-range tail as the terminal density is approached from below. 
By implication, for the special case of a thermal equilibrium process, the pair potential at 
the terminal density develops a long-range repulsive Coulombic form. The conclusion is that 
hyperuniformity at that terminal density is logically associated with the local electroneu- 
Tality condition that all equilibrium systems of electrostatically charged particles must obey 



The Ornstein-Zernike relation, though originally conceived to apply to systems in thermal 
equilibrium, can nevertheless be formally applied to any system for which the pair correlation 
function (72 (r) is available. Hyperuniform systems that are irregular and isotropic possess 
short-range pair correlation only, but as in the examples just cited the corresponding direct 
correlation functions are long-ranged. In an important sense, hyperuniform systems exhibit 
a kind of "inverted critical phenomenon". For conventional liquid- vapor critical points, 
h{r) = (72 (r) — 1 is long-ranged and implies diverging density fluctuations and isothermal 
compressibilities, while the direct correlation function c(r) remains short-ranged. Hyper- 
uniform systems have short range for h{r), vanishing volume fluctuations and isothermal 
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compressibility, and a long-ranged c(r). 

As a final matter, we mention that an attractive direction for future study of hyperuni- 
formity and related concepts involves consideration of collective density variables. These are 
defined by a nonlinear transformation of point-particle positions (1 < j < A^) as follows: 

N 

p(k) = 5:exp(zk-r,). (138) 

i=i 

If the particles interact through a spherically-symmetric pair potential whose Fourier trans- 
form exists and is denoted by V{k), then the overall potential energy for the N particles in 
volume V can be expressed in the following manner: 

* = ^E^(^)[p(k)p(-k)-iV]. (139) 

It has been demonstrated that at least in one dimension, application of a suitable V{k), 
followed by $ minimization, can totally suppress density fluctuations for k's near the origin. 
This automatically produces a hyperuniform system configuration. Analogous studies need 
to be pursued for two- and three-dimensional systems. 
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APPENDIX A: INTERSECTION VOLUME OF TWO IDENTICAL d- 
DIMENSIONAL SPHERES 

In this appendix, we obtain an explicit expression for the scaled intersection volume of 
two identical (i- dimensional spheres of radius R whose centers are separated by a distance 
r. This function a(r; R) is defined by (j2SI)- 



46 



We begin by noting that the d-dimensional Fourier transform ()16|) of any integrable 
function /(r) that depends only on the modulus r = |r| of the vector r is given by 2^ 

/»^,2.,f/;/-vw^Jf^<^. (Ai) 

and the inverse transform p7j] of f{k) is given by 



oo 



J(d/2)^i{kr) 



(27r)2 "'0 {^kry 

Here k is the modulus of the wave vector k and Ju{x) is the Bessel function of order z/. 
The Fourier transform of the window indicator function (jHT|) is given by 

/ 27r \ '^/^ 

= (^j R'JddkR). (A3) 
Therefore, the Fourier transform of a(r; R), defined by (f^ . is given by 

aik; R) = 2V/2r(l + d/2)t^^^^i^. (A4) 



Using the inverse transform ()A2|) yields the scaled intersection volume function to be 



r{d-2)/2 _/q J^d/2 

2 ' 2 



/i-.^f^,i)e(2i?-r), (A5) 



where 



is the normalized incomplete beta function (28l |. 

5,.(a,6)= r e-\l-tf-^dt, (A7) 

JO 



is the incomplete beta function, and 

r(a)r(&) 

/o " r(a + 6) 

is the beta function. 



B{a, h) = / - tf-^dt = \! (A8) 
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For the first five space dimensions, relation ()A5|1 . for r < 2R, yields 



a(r: R) = 1 , d 



a{r; R) = — 

TT 



a{r;R) = 1- 
2 

a{r] R) = — 



\2RJ 2R \ 



3 r 



.2 \ V2' 



TT 



COS 



a(r; i?) 
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16 ( 






r ^ 




' 5r 




-< 
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\ ^ 
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Figure ini shows graphs of the scaled intersection volume a(r; R) as a function of r for the first 
five space dimensions. For any dimension, a{r] R) is a monotonically decreasing function of 
r. At a fixed value of r in the open interval (0,2i?), a{r; R) is a monotonically decreasing 
function of the dimension d. 

1 



Spherical window of radius R 




0.4 0.8 1.2 1.6 

r/(2R) 



FIG. 9: The scaled intersection volume a{r;R) for spherical windows of radius i? as a function of 
r for the first five space dimensions. The uppermost curve is for d = 1 and lowermost curve is for 
d = 5. 

Expanding the general expression ()A5|) through first order in r for r < 2R yields 

where o{x) indicates terms of higher order than x. This relation will be of use to us in 
developing an asymptotic expression for the number variance for large windows. 
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APPENDIX B: FLUCTUATIONS IN EQUILIBRIUM HARD-PARTICLE SYS- 
TEMS 



Hard particles in equilibrium represent an example of a correlated system that is gen- 
erally not hyperuniform. The one-dimensional case of identical hard rods of length D in 
equilibrium is a particularly instructive case because the radial distribution function (72 (r) 
(in the thermodynamic limit) is known exactly for all densities js^: 



l2{x) = J2 Q{x-k) — 
J — 1 I r 



(x — k) 



k-l 



exp 



(x — k) 



(Bl) 



0)^(^-1)! 

where x = r/D is a dimensionless distance and = pD is the covering fraction of the rods, 
which lies in the closed interval [0, 1]. Below the close-packed space-filling value of = 1, the 
radial distribution function is a short-ranged function in the sense that one can always find 
a large enough value of r beyond which g2{r) remains appreciably close to unity. That is, 
for < 1, the correlation length is always finite. However, the point = 1 is singular in the 
sense that the system exhibits perfect long-range order and thus is hyperuniform. Indeed, 
at = 1, the nearest-neighbor distance for each rod is exactly equal to D: a situation that 
is identically the same as the single-scale one-dimensional periodic point pattern studied in 
Section IV. 

1 



Equilibrium Hard Rods 




FIG. 10: The volume coefficient ^(0) = (1 — 0)'^ and surface-area coefficient -6(0) [defined in 
relation (|35() ] as a function of the reduced density for a one-dimensional system of identical hard 
rods in equilibrium. At the hyperuniform density (j) = 1, B = A/2 = 1/12, which corresponds to 
the perfectly ordered close-packed state. 



49 



Using relation ()B1|) in conjunction with relations ()36p and ()37|) enables us to compute the 
"volume" and "surface-area" contributions to the variance as a function of reduced density 
for identical hard rods in equilibrium. The results are summarized in Fig. ^1 We see that 
as the density increases, the volume fluctuations decrease monotonically and only vanish at 
the space-filling density 0=1: the hyperuniform state. Of course, B vanishes at = 
and increases in value as increases until it achieves a maximum value at ~ 0.5. At the 
hyperuniform state (0 = 1), B = A/2 = 1/12, which corresponds to the perfectly ordered 
close-packed state. For sufficiently small densities, the surface-area coefficient of equilibrium 
hard-sphere systems in higher dimensions is expected to have the same qualitative behavior 
as the one-dimensional case. Specifically, the same trends should occur in higher dimensions 
for densities in the range < < 0/, where 0/ corresponds to the freezing density, i.e., the 
point above which the system undergoes a disorder to order phase transition. For densities 
between freezing and melting points, the behavior of the surface-area coefficient is expected 
to be qualitatively different from that for hard rods in equilibrium, which is devoid of a phase 
transition. However, we can definitively assert that the highest achievable density along the 
stable crystal branch is a hyperuniform state. In particular, for hard disks {d = 2) and 
hard spheres {d = 3) in equilibrium, the hyperuniform states correspond to the close-packed 
triangular lattice and FCC lattice, respectively. 

APPENDIX C: HOW SMALL CAN THE VOLUME COEFFICIENT BE FOR 
HYPOSURFICIAL SYSTEMS? 

We know that a statistically homogeneous and isotropic point pattern cannot simul- 
taneously be hyperuniform and hyposurficial, i.e., the volume coefficient A [cf. (I43|l ] and 
surface-area coefficient B [cf. fl44j) ] cannot both be zero for a strictly convex window (Section 
IIC). The purpose of this appendix is to investigate how small A can be made for an infinite 
hyposurficial point pattern {B = 0). To that end we consider a hypothetical spherically 
symmetric pair correlation function (72 (r) and a spherical window. We do not place any 
additional restrictions on g2{r) besides the necessary realizability conditions that g2{r) > 
for all r and S{k) > for all k. The hypothetical correlation function is characterized by 
three parameters e, C, and D as follows: 

92ir) = gs{r)+gL{r), (CI) 
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where gsi^) denotes the short-ranged part defined by the step function 

<r < D, 

9s{r) = { - - - ' (C2) 

r > D, 

and giir) denotes the long-ranged part defined by 

0, 0<r<L', 

9L{r) = { /Dx'^+i+^ (C3) 
' Ce - , r>D, 



r ) 

Here D is a length parameter, C is a dimensionless constant, and e is a positive (e > 0) but 
small parameter. The necessary condition g2{T) > requires that the constant C satisfy the 
trivial inequality 

C > -e~\ (C4) 

The form of g2 ensures that we can make the surface-area coefficient B vanish identically, 
as required. According to relation (0H), the surface-area coefficient B is proportional to the 
dth moment of the total correlation function h{r) = g2{r) — 1. The dth moment integral for 
the hypothetical pair correlation function I^Cl^ is given by 



hhydr = + CD'^+\ (C5) 

^ ^ d + 1 



To make this integral vanish, we take 



C=jl^>0, (C6) 



which of course satisfies the inequality (jCMjl . For such a hyposurficial correlation function 
that also satisfies the nonnegativity condition S{k) > 0, we now show that the volume 
coefficient A is only nonzero by 0{e^). 

Consider volume coefficient A [cf. with this value of C: 

A= \im S{k) = l + p [ h{r)dr = l-2''<j)+(^]-^, (C7) 
|k|^o ^7sRd ^ \d+ljl + e ' 

where = pvi{D/2) is a dimensionless density. If one incorrectly sets A to be zero, one 
finds that the corresponding density is given by 

1 



2'^ I 1 



de 



(C8) 



(rf+l)(l + e)_ 
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At such a value of 0, however, S{k) will be negative for some k > near the origin for 
sufficiently small but nonzero e, which shows in this specific instance that the point pat- 
tern corresponding to such a hypothetical g2 cannot simultaneously be hyperuniform and 
hyposurficial, as expected. However, one can make S{k = 0) positive and very small (while 
satisfying S{k) > for all k) at a value slightly smaller than ()C8|) in the limit e — 0+. 

The other necessary condition S{k) > will be obeyed for all k provided that the number 
density is no larger than some "terminal density" pc (or (pc) jsQ, 40, 41 1. The structure factor 
is given by 

S{k) = l+p[Hs{k) + HL{k)], (C9) 

where Hs{k) and Hi{k) are the Fourier transforms of gsi^) — 1 and gL{f)i respectively. The 
terminal density is given by 

= ~miMHs{k) + HL{k)\ 

For simplicity, we will specialize to the case d = 3, keeping in mind that our general 
conclusions apply to arbitrary dimension. Based on the aforementioned arguments, it is 
sufficient to consider the behavior of S{k) for small k: 



S{k) 



-1 + ^ - 0\{kDY] 
30v/^e r(i-f),,^^„^^ 



+ 



1 + e 



kD\'+' + -^{kDY + 0[{kDY]. (Cll) 



21+^(1 + e)r(2 + f)' ' 1-e 

The nonanalytic term \kD\^'^^ [which arises due to inclusion of relation ()C3|) for gi] has the 
effect of displacing the minimum of S{k) away from the origin when (7^ = to a symmetric 
pair of locations determined by 

\km^nD\ = (C12) 
lo 

as e — s> 0^. Moreover, in this leading order 

5(0) - S{krmn) = (C13) 

o4 

Note that this would lead to an 0{e^) correction to expression (jCvSIl for 0*. In summary, by 
adopting the correlation function ^CA\ with C = l/{d + 1), we can make the surface-area 
coefficient B = and at the terminal density (pc, the structure factor S{0) = A is only 
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nonzero by 0{e^). 
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